#pragma once

#include "vector.h"
#include <algorithm>
#include <functional>
//for IO operations #include<glib/g3/io.h>

namespace g3 {

template <class T>
class basic_matrix;

typedef basic_matrix<double> matrix;
typedef basic_matrix<float> matrixf;
typedef basic_matrix<long double> matrixl;

template <class T>
class basic_matrix {
public:
    T e[3][3];
public:
	typedef T * iterator;
	typedef const T * const_iterator;
public:
	T * begin()
		{ return &e[0][0]; }
	const T * begin() const
		{ return &e[0][0]; }
	T * end()
		{ return &e[2][2] + 1; }
	const T * end() const
		{ return &e[2][2] + 1; }
	static int size() 
		{ return 9; }
public:
	basic_matrix()
		{ }
	template <class U> 
	basic_matrix(const basic_matrix<U> & v) 
		{ std::copy(v.begin(), v.end(), begin()); }
	basic_matrix & identity(T diagval = 1) {
		e[0][0] = e[1][1] = e[2][2] = diagval;
		e[0][1] = e[0][2] = e[1][0] = e[1][2] = e[2][0] = e[2][1] = 0;
		return * this;
	}
	basic_matrix & zero() 
		{ std::fill(begin(), end(), T(0)); return * this; } 
public:
	const T & operator ()(int i, int j) const
		{ return e[i][j]; }
	T & operator ()(int i, int j)
		{ return e[i][j]; }
	/// multiplication product of the matrix (lhs) and given vector (rhs)
	basic_vector<T> operator ()(const basic_vector<T> & v) const {
		return basic_vector<T>(scalar_product(row(0), v),
			scalar_product(row(1), v), scalar_product(row(2), v));
	}
	/// multiplication product of given vector (lhs) and the matrix (rhs)
	basic_vector<T> left_product(const basic_vector<T> & v) const {
		return basic_vector<T>(scalar_product(col(0), v),
			scalar_product(col(1), v), scalar_product(col(2), v));
	}
public:
	basic_vector<T> col(int i) const
		{ return basic_vector<T>(e[0][i], e[1][i], e[2][i]); }
	basic_matrix<T> & set_col(int i, const basic_vector<T> & v)
		{ e[0][i] = v.x; e[1][i] = v.y; e[2][i] = v.z; return * this; }
	basic_vector<T> row(int i) const
		{ return basic_vector<T>(e[i][0], e[i][1], e[i][2]); }
	basic_matrix<T> & set_row(int i, const basic_vector<T> & v)
		{ e[i][0] = v.x; e[i][1] = v.y; e[i][2] = v.z; return * this; }
	T det() const
		{ return mixed_product(col(0), col(1), col(2)); };
	operator bool() const ///< true if not degenerate
		{ return det() != 0; };
	basic_matrix<T> & transpose(); ///< transposes itself
	basic_matrix<T> transposed() const; ///< returns transposed matrix, while leaves itself intact
	basic_matrix<T> inverse() const;
	basic_vector<T> inverse(const basic_vector<T> & v) const; ///< \return product of inverse matrix and given vector
	T tr() const ///< trace of the matrix
		{ return e[0][0] + e[1][1] + e[2][2]; }
public:
    basic_matrix & operator += (const basic_matrix<T> & v)
		{ std::transform(begin(), end(), v.begin(), begin(), std::plus<T>()); return * this; }
    basic_matrix & operator += (T a) ///< adds identity matrix multiplied on a
		{ e[0][0] += a; e[1][1] += a; e[2][2] += a; return * this; }
    basic_matrix & operator -= (const basic_matrix<T> & v)
		{ std::transform(begin(), end(), v.begin(), begin(), std::minus<T>()); return * this; }
    basic_matrix & operator -= (T a) ///< adds identity matrix multiplied on a
		{ e[0][0] -= a; e[1][1] -= a; e[2][2] -= a; return * this; }
    basic_matrix & operator *= (T a)
		{ std::transform(begin(), end(), begin(), std::bind2nd(std::multiplies<T>(), a)); return * this; }
    basic_matrix & operator /= (T a)
		{ std::transform(begin(), end(), begin(), std::bind2nd(std::divides<T>(), a)); return * this; }
	bool operator == (const basic_matrix<T> & v) const
		{ return std::equal(begin(), end(), v.begin()); }
	bool operator != (const basic_matrix<T> & v) const
		{ return !(*this == v); }
};

//inlines

template <class T>
inline basic_matrix<T> & basic_matrix<T>::transpose()
{
	std::swap(e[0][1], e[1][0]);
	std::swap(e[0][2], e[2][0]);
	std::swap(e[1][2], e[2][1]);
	return * this;
}

template <class T>
inline basic_matrix<T> basic_matrix<T>::transposed() const
{
	return basic_matrix<T>(*this).transpose();
}

template <class T>
inline basic_matrix<T> basic_matrix<T>::inverse() const {
	return basic_matrix<T>()
		.set_row(0, outer_product(col(1), col(2)))
		.set_row(1, outer_product(col(2), col(0)))
		.set_row(2, outer_product(col(0), col(1))) / det();
}

template <class T>
inline basic_vector<T> basic_matrix<T>::inverse(const basic_vector<T> & v) const {
	return basic_vector<T>(
		mixed_product(v, col(1), col(2)),
		mixed_product(v, col(2), col(0)),
		mixed_product(v, col(0), col(1))
	) / det();
}

template <class T>
inline basic_matrix<T> operator + (const basic_matrix<T> & a, const basic_matrix<T> & b) {
	basic_matrix<T> t(a);
    return t += b;
}

template <class T>
inline basic_matrix<T> operator + (const basic_matrix<T> & v, T a) {
	basic_matrix<T> t(v);
    return t += a;
}

template <class T>
inline basic_matrix<T> operator - (const basic_matrix<T> & a, const basic_matrix<T> & b) {
	basic_matrix<T> t(a);
    return t -= b;
}

template <class T>
inline basic_matrix<T> operator - (const basic_matrix<T> & v, T a) {
	basic_matrix<T> t(v);
    return t -= a;
}

template <class T>
inline basic_matrix<T> operator * (const basic_matrix<T> & v, T a) {
	basic_matrix<T> t(v);
    return t *= a;
}

template <class T>
inline basic_matrix<T> operator * (T a, const basic_matrix<T> & v) {
	basic_matrix<T> t(v);
    return t *= a;
}

template <class T>
inline basic_vector<T> operator * (const basic_matrix<T> & m, const basic_vector<T> & v) {
    return m(v);
}

template <class T>
inline basic_vector<T> operator * (const basic_vector<T> & v, const basic_matrix<T> & m) {
    return m.left_product(v);
}

template <class T>
inline basic_matrix<T> operator * (const basic_matrix<T> & a, const basic_matrix<T> & b) {
    basic_matrix<T> t;
	for (int i = 0; i < 3; ++i)
		for (int j = 0; j < 3; ++j)
			t(i,j) = scalar_product(a.row(i), b.col(j));
	return t;
}

template <class T>
inline basic_matrix<T> operator / (const basic_matrix<T> & v, T a) {
	basic_matrix<T> t(v);
    return t /= a;
}

} //namespace g3
